Variable type

First, a descriptive exploration of the variables was performed in order to understand them better and detect possible problems.

skim(data)

Variable type: factor

Variable type: numeric

After visually analyzing the distribuitions of all the features, no abnormalities seem to be found. Other than the target variable, which is a factor, all variables are numeric. The type of variable is correctly assigned to all features.
In addition, there are no missing values.

Dataset features

ggplot(data, aes(diagnosis, fill=diagnosis)) + 
  geom_bar() +
  labs(x="Diagnosis", y="Number of patients") +
  guides(fill=FALSE) +
  scale_fill_manual( values = c("#686868","#9F2042"))

library(gridExtra)
p <- list()
for (j in colnames(data)[2:11]) {
  p[[j]] <- ggplot(data=data, aes_string(x="diagnosis", y=j)) + 
            geom_boxplot(aes(fill=factor(diagnosis))) + guides(fill=FALSE) +
            theme(axis.title.y = element_text(size=8)) +  
            geom_jitter(alpha = 0.2, width = 0.2) +
            scale_fill_manual( values = c("#686868","#9F2042"))
}
do.call(grid.arrange, c(p, ncol=5))

p <- list()
for (j in colnames(data)[12:21]) {
  p[[j]] <- ggplot(data=data, aes_string(x="diagnosis", y=j)) + 
            geom_boxplot(aes(fill=factor(diagnosis))) + guides(fill=FALSE) +
            theme(axis.title.y = element_text(size=8)) +  
            geom_jitter(alpha = 0.2, width = 0.2) +
            scale_fill_manual( values = c("#686868","#9F2042"))
}
do.call(grid.arrange, c(p, ncol=5))

p <- list()
for (j in colnames(data)[22:31]) {
  p[[j]] <- ggplot(data=data, aes_string(x="diagnosis", y=j)) + 
            geom_boxplot(aes(fill=factor(diagnosis))) + guides(fill=FALSE) +
            theme(axis.title.y = element_text(size=8)) +  
            geom_jitter(alpha = 0.2, width = 0.2) +
            scale_fill_manual( values = c("#686868","#9F2042"))
}
do.call(grid.arrange, c(p, ncol=5))

library(pheatmap)
pheatmap(cor(data[,-1]))

Stratified sampling

As seen before, there are more benign that malignan samples. Stratified sampling was used in order to guarantee the same proportion in each class than the one had in the complete data set.

library(rsample)
set.seed(123)
split_strat <- initial_split(data,prop = 0.8, strata = "diagnosis")

datos_train <- training(split_strat)
datos_test <- testing(split_strat)

prop.table(table(data$diagnosis))

## 
##         B         M 
## 0.6274165 0.3725835

table(datos_test$diagnosis) %>% prop.table()

## 
##         B         M 
## 0.6283186 0.3716814

table(datos_train$diagnosis) %>% prop.table()

## 
##        B        M 
## 0.627193 0.372807

Imputation of missing values

There are no missing values, thus there is no need to imputate.

Variables with variance close to zero

As seen on the table below, all predictors have a significant variance. Hence, no features were removed based on this criteria.

datatable(nearZeroVar(data,saveMetrics = T), class = 'cell-border stripe') 

Standardization and scaling

The boxplots demonstrated the need for centering and scaling. The difference in scales can have a great impact in the model. The step_normalize function centers and scales the data.

obj_recipe <- obj_recipe %>%
  step_normalize(all_numeric()) 

obj_recipe

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         30
## 
## Operations:
## 
## Centering and scaling for all_numeric()

Once the recipe object has been created with the preprocessing transformations, they are learned with the training data and applied to the two sets.

trained_recipe <- obj_recipe %>%  prep(training=datos_train)
trained_recipe

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         30
## 
## Training data contained 456 data points and no missing data.
## 
## Operations:
## 
## Centering and scaling for radius_mean, texture_mean, ... [trained]

datos_train_prep <-trained_recipe %>% bake(new_data = datos_train)
datos_test_prep <- trained_recipe %>% bake(new_data = datos_test)

A “glimpse” of the values of the features after applying the transformations is shown below.

glimpse(datos_train_prep)

## Rows: 456
## Columns: 31
## $ radius_mean             <dbl> 1.1026662, 1.8353345, 1.5854321, -0.7630822, …
## $ texture_mean            <dbl> -2.0368159, -0.3414951, 0.4568427, 0.2572583,…
## $ perimeter_mean          <dbl> 1.27735149, 1.69374265, 1.57418479, -0.586926…
## $ area_mean               <dbl> 0.9838048, 1.9037595, 1.5555920, -0.7567497, …
## $ smoothness_mean         <dbl> 1.59709014, -0.82451799, 0.96398998, 3.330921…
## $ compactness_mean        <dbl> 3.33607335, -0.48595341, 1.07505329, 3.457096…
## $ concavity_mean          <dbl> 2.69220042, -0.01326369, 1.38895856, 1.947309…
## $ concave.points_mean     <dbl> 2.58556012, 0.57146330, 2.08288679, 1.4885803…
## $ symmetry_mean           <dbl> 2.30706285, 0.02186195, 0.98940170, 2.9771876…
## $ fractal_dimension_mean  <dbl> 2.31733476, -0.88939910, -0.40635208, 5.04247…
## $ radius_se               <dbl> 2.37985886, 0.47285546, 1.17168663, 0.3072245…
## $ texture_se              <dbl> -0.52600515, -0.83843605, -0.74182673, -0.069…
## $ perimeter_se            <dbl> 2.72216499, 0.25563868, 0.81964691, 0.2779709…
## $ area_se                 <dbl> 2.325575265, 0.686914706, 1.099058909, -0.280…
## $ smoothness_se           <dbl> -0.19032973, -0.58560573, -0.27416595, 0.7224…
## $ compactness_se          <dbl> 1.356680925, -0.683695689, 0.847154173, 2.805…
## $ concavity_se            <dbl> 0.71209089, -0.41400693, 0.21812085, 0.804409…
## $ concave.points_se       <dbl> 0.6825459, 0.2843963, 1.4417704, 1.1338896, -…
## $ symmetry_se             <dbl> 1.172419370, -0.779156596, 0.261925750, 4.751…
## $ fractal_dimension_se    <dbl> 0.913959900, -0.083344862, 0.306057486, 2.043…
## $ radius_worst            <dbl> 1.8696583, 1.7894970, 1.4976275, -0.2823654, …
## $ texture_worst           <dbl> -1.34090113, -0.35703897, -0.01398177, 0.1429…
## $ perimeter_worst         <dbl> 2.30085932, 1.53384537, 1.34655126, -0.247826…
## $ area_worst              <dbl> 1.96064959, 1.85181528, 1.42511569, -0.546512…
## $ smoothness_worst        <dbl> 1.338856915, -0.379130952, 0.542497956, 3.468…
## $ compactness_worst       <dbl> 2.69410378, -0.43561800, 1.11878871, 4.005450…
## $ concavity_worst         <dbl> 2.15333090, -0.14259294, 0.87673281, 2.031285…
## $ concave.points_worst    <dbl> 2.33628782, 1.11076543, 1.99054851, 2.2143529…
## $ symmetry_worst          <dbl> 2.8082563, -0.2375773, 1.1824953, 6.1601542, …
## $ fractal_dimension_worst <dbl> 2.02437558, 0.29713613, 0.21389567, 5.1516732…
## $ diagnosis               <fct> M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, …

p <- list()
for (j in colnames(datos_train_prep)[1:10]) {
  p[[j]] <- ggplot(data=datos_train_prep, aes_string(x="diagnosis", y=j)) + 
            geom_boxplot(aes(fill=factor(diagnosis))) + guides(fill=FALSE) +
            theme(axis.title.y = element_text(size=8)) +  
            geom_jitter(alpha = 0.2, width = 0.2) +
            scale_fill_manual( values = c("#686868","#9F2042"))
}
do.call(grid.arrange, c(p, ncol=5))

Model 1: Glmnet

Glmnet fits a generalized linear model via penalized maximum likelihood. The penalization parameters are alpha and lambda, the hypeparameter selection was left by default. These hyperparameters deal with correlated predictors, which is of importance in this project.

set.seed(123)
modelo_glm <- train(diagnosis ~.,
                    method="glmnet", 
                    family="binomial", 
                    trControl=control_train, 
                    data=datos_train_prep, 
                    metric="Accuracy")
modelo_glm$bestTune

##   alpha      lambda
## 5  0.55 0.007608795

ggplot(modelo_glm, highlight = T)

The most important variable is radius_worst.

vip(modelo_glm, geom = "point")

confMat(modelo_glm)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 71  0
##          M  0 42
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9679, 1)
##     No Information Rate : 0.6283     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##                                      
##  Mcnemar's Test P-Value : NA         
##                                      
##             Sensitivity : 1.0000     
##             Specificity : 1.0000     
##          Pos Pred Value : 1.0000     
##          Neg Pred Value : 1.0000     
##              Prevalence : 0.6283     
##          Detection Rate : 0.6283     
##    Detection Prevalence : 0.6283     
##       Balanced Accuracy : 1.0000     
##                                      
##        'Positive' Class : B          
## 

Model 2: Random Forest

The Random forest is based on constructing a multitude of decision trees.

set.seed(123)
hip <- data.frame(mtry=1:30) # Randomly selected predictors
modelo_rf <-train(diagnosis ~.,
                  method="rf", 
                  trControl=control_train, 
                  data=datos_train_prep, 
                  tuneGrid=hip,  
                  metric="Accuracy")

modelo_rf$bestTune

##   mtry
## 1    1

ggplot(modelo_rf, highlight = T)

The most important variable of the model is area_worst.

vip(modelo_rf, geom = "point")

confMat(modelo_rf)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 69  3
##          M  2 39
##                                           
##                Accuracy : 0.9558          
##                  95% CI : (0.8998, 0.9855)
##     No Information Rate : 0.6283          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9048          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.9718          
##             Specificity : 0.9286          
##          Pos Pred Value : 0.9583          
##          Neg Pred Value : 0.9512          
##              Prevalence : 0.6283          
##          Detection Rate : 0.6106          
##    Detection Prevalence : 0.6372          
##       Balanced Accuracy : 0.9502          
##                                           
##        'Positive' Class : B               
## 

Model 3: AdaBoost

The next algorithm used was AdaBoost. The hypeparameter selection for AdaBoost was left by default, one of the two methods and the number of iterations can be chosen.

set.seed(123)
modelo_ada <-train(diagnosis ~.,
                   method="adaboost", 
                   trControl=control_train, 
                   data=datos_train_prep, 
                   metric="Accuracy")
modelo_ada$bestTune

##   nIter        method
## 2    50 Real adaboost

ggplot(modelo_ada, highlight = T)

confMat(modelo_ada)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 70  9
##          M  1 33
##                                           
##                Accuracy : 0.9115          
##                  95% CI : (0.8433, 0.9567)
##     No Information Rate : 0.6283          
##     P-Value [Acc > NIR] : 6.062e-12       
##                                           
##                   Kappa : 0.8029          
##                                           
##  Mcnemar's Test P-Value : 0.02686         
##                                           
##             Sensitivity : 0.9859          
##             Specificity : 0.7857          
##          Pos Pred Value : 0.8861          
##          Neg Pred Value : 0.9706          
##              Prevalence : 0.6283          
##          Detection Rate : 0.6195          
##    Detection Prevalence : 0.6991          
##       Balanced Accuracy : 0.8858          
##                                           
##        'Positive' Class : B               
## 

Model 4: Support Vector Machines with Linear and Polynomial Kernel

SVM training algorithm builds a model that assigns new examples to one category or the other. The classification performed can be linear or non-linear. In this case, the hyperparameters are cost, degree and scale. The chosen values are shown on the chunk bellow. The best model was the linear model.

set.seed(123)
hip_svmP <- expand.grid(C=c(0.001, 0.01, 0.1, 0.5, 1, 10),degree=c(1,2,3),scale=1)
modelo_svmPoly <- train(diagnosis ~.,
                    method = "svmPoly",
                    trControl = control_train,
                    data = datos_train_prep,
                    tuneGrid = hip_svmP,
                    metric = "Accuracy")
modelo_svmPoly$bestTune

##   degree scale   C
## 7      1     1 0.1

ggplot(modelo_svmPoly, highlight = T)

confMat(modelo_svmPoly)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 71  0
##          M  0 42
##                                      
##                Accuracy : 1          
##                  95% CI : (0.9679, 1)
##     No Information Rate : 0.6283     
##     P-Value [Acc > NIR] : < 2.2e-16  
##                                      
##                   Kappa : 1          
##                                      
##  Mcnemar's Test P-Value : NA         
##                                      
##             Sensitivity : 1.0000     
##             Specificity : 1.0000     
##          Pos Pred Value : 1.0000     
##          Neg Pred Value : 1.0000     
##              Prevalence : 0.6283     
##          Detection Rate : 0.6283     
##    Detection Prevalence : 0.6283     
##       Balanced Accuracy : 1.0000     
##                                      
##        'Positive' Class : B          
## 

Model 5: Naive Bayes

The Naive Bayes algorithm is based on applying Bayes’ theorem with strong (naïve) independence assumptions between the features. The hypermarameters were left by default. The results obtained are not as favorable compared to the other models probably because it doesn’t consider any kind of correlation between features.

set.seed(123)
modelo_nb <- train(diagnosis ~.,
                    method = "nb",
                    trControl = control_train,
                    data = datos_train_prep,
                    metric = "Accuracy")
modelo_nb$bestTune

##   fL usekernel adjust
## 2  0      TRUE      1

ggplot(modelo_nb, highlight = T)

confMat(modelo_nb)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 66  2
##          M  5 40
##                                           
##                Accuracy : 0.9381          
##                  95% CI : (0.8765, 0.9747)
##     No Information Rate : 0.6283          
##     P-Value [Acc > NIR] : 1.718e-14       
##                                           
##                   Kappa : 0.8693          
##                                           
##  Mcnemar's Test P-Value : 0.4497          
##                                           
##             Sensitivity : 0.9296          
##             Specificity : 0.9524          
##          Pos Pred Value : 0.9706          
##          Neg Pred Value : 0.8889          
##              Prevalence : 0.6283          
##          Detection Rate : 0.5841          
##    Detection Prevalence : 0.6018          
##       Balanced Accuracy : 0.9410          
##                                           
##        'Positive' Class : B               
## 

Model 6: k-Nearest Neighbors

K-nearest neighbors is an algorithm that classifies new cases based on a similarity measure (distance between neighbors). The hyperparameters were also left by default.

set.seed(123)
modelo_kknn <- train(diagnosis ~.,
                    method = "kknn",
                    trControl = control_train,
                    data = datos_train_prep,
                    metric = "Accuracy")
modelo_kknn$bestTune

##   kmax distance  kernel
## 2    7        2 optimal

ggplot(modelo_kknn, highlight = T)

confMat(modelo_kknn)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  B  M
##          B 70  1
##          M  1 41
##                                           
##                Accuracy : 0.9823          
##                  95% CI : (0.9375, 0.9978)
##     No Information Rate : 0.6283          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.9621          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.9859          
##             Specificity : 0.9762          
##          Pos Pred Value : 0.9859          
##          Neg Pred Value : 0.9762          
##              Prevalence : 0.6283          
##          Detection Rate : 0.6195          
##    Detection Prevalence : 0.6283          
##       Balanced Accuracy : 0.9811          
##                                           
##        'Positive' Class : B               
## 

Comparison of first models

The accuracy and kappa metrics of the 10-folds for each model is shown on the following table.

modelos <- list(GLM=modelo_glm, RF=modelo_rf, ADA=modelo_ada, SVM=modelo_svmPoly, NB=modelo_nb, KNN=modelo_kknn)
results_resamples <- resamples(modelos)

datatable(results_resamples$values)

The time taken to compute all the folds and the final model for each algorithm is shown on the next table.

datatable(results_resamples$timings)

Some models obtained an accuracy of 1 on one of the folds. For this reason, in order to decide which model had a better performance the mean accuracy of all folds was computed.

metricas_resamples <- results_resamples$values %>%
                         gather(key = "modelo", value = "valor", -Resample) %>%
                         separate(col = "modelo", into = c("modelo", "metrica"), sep = "~", remove = TRUE)

metricas_resamples <- metricas_resamples %>% filter(metrica == "Accuracy") %>%
  group_by(modelo) %>% 
  mutate(media = mean(valor)) %>%
  ungroup() %>%
  ggplot(aes(x = reorder(modelo, media), y = valor, color = modelo)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.1, alpha = 0.6) +
    theme_bw() +
    labs(title = "Validation: Mean Accuracy of the repeated-CV",
         subtitle = "Models are ordered based on the mean") +
    coord_flip()

metricas_resamples

The models with the highest accuracy are the SVM model and the ADA model. The SVM algorithm has the highest mean and takes less time to execute than the ADA algorithm, which is actually the one that takes the longest. Their ROC curves are shown below, the AUC of both models demonstrate their perfect performance.

library(ROCR)
predictions <- predict(modelo_ada, datos_test_prep %>% dplyr::select(-diagnosis))
pred <- prediction(as.numeric(predictions),as.numeric(datos_test_prep$diagnosis))
perf <- performance(pred,"tpr","fpr")
AUC <- as.numeric(performance(pred,"auc")@y.values)

plot(perf,colorize=TRUE)
text(0.5, 0.4, paste("AUC = ",round(AUC,4)))
title("ADA")
grid()

predictions <- predict(modelo_svmPoly, datos_test_prep %>% dplyr::select(-diagnosis))
pred <- prediction(as.numeric(predictions),as.numeric(datos_test_prep$diagnosis))
perf <- performance(pred,"tpr","fpr")
AUC <- as.numeric(performance(pred,"auc")@y.values)

plot(perf,colorize=TRUE)
text(0.5, 0.4, paste("AUC = ",round(AUC,4)))
title("SVM")
grid()

Comparison of models with PCA

library(rlist)
modelos2 <- list.append(modelos, GLMpca=modelo_glm2, RFpca=modelo_rf2, ADApca=modelo_ada2, SVMpca=modelo_svmPoly2, NBpca=modelo_nb2, KNNpca=modelo_kknn2)
results_resamples2 <- resamples(modelos2)

Metrics…

datatable(results_resamples2$values)

Time…

datatable(results_resamples2$timings)

metricas_resamples2 <- results_resamples2$values %>%
                         gather(key = "modelo", value = "valor", -Resample) %>%
                         separate(col = "modelo", into = c("modelo", "metrica"), sep = "~", remove = TRUE)

metricas_resamples2 <- metricas_resamples2 %>% filter(metrica == "Accuracy") %>%
  group_by(modelo) %>% 
  mutate(media = mean(valor)) %>%
  ungroup() %>%
  ggplot(aes(x = reorder(modelo, media), y = valor, color = modelo)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.1, alpha = 0.6) +
    theme_bw() +
    labs(title = "Validation: Mean Accuracy of the repeated-CV",
         subtitle = "Models are ordered based on the mean") +
    coord_flip()

As seen on the graph below, the SVM model and the ADA model still outperform all other models. The first SVM model still has the highest accuracy mean indicating that reducing dimensionality does not necessarily improve the model, regardless of the number of variables or the correlation. If execution time is taken into account the GLMpca model is much faster and also has a good performance.

metricas_resamples2

Finally, the top two algorithms (SVM and ADA) were again used to create two models that take into account the class imbalance. All the features were taken into account.

Final comparison

modelos3 <- list.append(modelos2, ADAweights=modelo_ada3,SVMweights=modelo_svmPoly3)
results_resamples3 <- resamples(modelos3)

Metrics…

datatable(results_resamples3$values)

Time…

datatable(results_resamples3$timings)

metricas_resamples3 <- results_resamples3$values %>%
                         gather(key = "modelo", value = "valor", -Resample) %>%
                         separate(col = "modelo", into = c("modelo", "metrica"), sep = "~", remove = TRUE)

metricas_resamples3 <- metricas_resamples3 %>% filter(metrica == "Accuracy") %>%
  group_by(modelo) %>% 
  mutate(media = mean(valor)) %>%
  ungroup() %>%
  ggplot(aes(x = reorder(modelo, media), y = valor, color = modelo)) +
    geom_boxplot(alpha = 0.6, outlier.shape = NA) +
    geom_jitter(width = 0.1, alpha = 0.6) +
    theme_bw() +
    labs(title = "Validation: Mean Accuracy of the repeated-CV",
         subtitle = "Models are ordered based on the mean") +
    coord_flip()

metricas_resamples3

Surprisingly, the SVM model that takes into account the proportion of samples in each class has the same accuracy mean as the SVM model but takes a bit less time to compute.